library(Seurat)
library(SingleCellExperiment)
library(DropletQC)
library(scDblFinder)
library(ggplot2)
library(qpdf)
library(openxlsx)

proj_path <- dirname(getwd())

WT_dir <- file.path(getwd(), "WT_QC")
dir.create(WT_dir)

WT_plots <- file.path(WT_dir, "WT_Plots")
dir.create(WT_plots)

WT_rds <- file.path(WT_dir, "WT_RDS_Files")
dir.create(WT_rds)

WT_markers <- file.path(WT_dir, "WT_Markers")
dir.create(WT_markers)
########################### Create WT Seurat Object ############################

# Import count matrix for the wildtype 
WT_mat <- Read10X(file.path(proj_path, 
                            "WT",
                            "outs",
                            "filtered_feature_bc_matrix"),  
                  gene.column = 2,  
                  cell.column = 1,  
                  unique.features = TRUE,  
                  strip.suffix = FALSE)

# Create seurat object for the wildtype
WT_s <- CreateSeuratObject(WT_mat, 
                           project = "WT")

# Calculate percentages of mitochondrial and ribosomal genes
WT_s <- PercentageFeatureSet(WT_s, 
                             "^mt-", 
                             col.name = "pct_mito")
WT_s <- PercentageFeatureSet(WT_s, 
                             "^Rp[sl]\\d+", 
                             col.name = "pct_ribo")

# Normalize 
WT_s <- NormalizeData(WT_s)

# Number of cells
nrow(WT_s@meta.data) # 15953
################################# WT DropletQC #################################

WT_nf <- nuclear_fraction_tags(outs = file.path(proj_path,
                                                "WT",
                                                "outs"),
                               tiles = 1, cores = 9, verbose = FALSE)

saveRDS(WT_nf, file.path(WT_rds, "0_WT_nuclear_fraction.rds"))

identical(row.names(WT_s@meta.data), row.names(WT_nf)) # True

WT_s <- AddMetaData(WT_s, WT_nf$nuclear_fraction, col.name = "nuclear_fraction")

################################################################################

WT_s.nf.umi <- data.frame(nf=WT_s$nuclear_fraction,
                         umi=WT_s$nCount_RNA)

# Identify empty droplets
WT_s.ed <- identify_empty_drops(nf_umi=WT_s.nf.umi)

################################################################################

WT_s <- FindVariableFeatures(WT_s, 
                             selection.method = "vst", 
                             nfeatures = 2000)
WT_s <- ScaleData(WT_s)
WT_s <- RunPCA(WT_s)

# Determine number of dimensions that capture at least 80% of the variation
pct <- WT_s[["pca"]]@stdev / sum(WT_s[["pca"]]@stdev) * 100
cumu <- cumsum(pct)
PC_80_pct <- which(cumu >= 80)[1]
PC_80_pct #32

# UMAP and clustering
WT_s <- FindNeighbors(WT_s, 
                      dims = 1:PC_80_pct)
WT_s <- RunUMAP(WT_s, 
                dims=1:PC_80_pct)

WT_s <- FindClusters(WT_s, 
                     resolution = 0.6)

WT_s <- AddMetaData(WT_s, WT_s$seurat_clusters, col.name = "WT_SC_1")

################################################################################

WT_s.ed$cell_type <- WT_s$WT_SC_1

# Identify damaged cells
WT_s.ed.dc <- identify_damaged_cells(WT_s.ed, verbose = FALSE, output_plots = FALSE)

WT_s <- AddMetaData(WT_s, WT_s.ed.dc$df$cell_status, col.name = "DropletQC")

# Save the seurat object
saveRDS(WT_s, file.path(WT_rds, "0_WT_seurat.rds"))
############################# Unfiltered WT Plots ##############################

# Ranges of counts, expressed genes, and percent mitochondria gene expression
range(WT_s$nCount_RNA) # min = 500; max = 101497
range(WT_s$nFeature_RNA) # min = 36; max = 10055
range(WT_s$pct_mito) # min = 0.04%; max = 96.42%
range(WT_s$pct_ribo) # min = 0.11%; max = 70.49%


pdf(file.path(WT_plots, "0_WT_unfiltered_scatter_plots.pdf"), height = 8, width = 12)
FeatureScatter(WT_s,
               feature1 = "nCount_RNA", 
               feature2 = "nFeature_RNA", 
               group.by = "DropletQC") +
  ggtitle("WT - Unfiltered")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_s,
               feature1 = "nCount_RNA", 
               feature2 = "pct_mito", 
               group.by = "DropletQC") +
  ggtitle("WT - Unfiltered")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_s,
               feature1 = "nCount_RNA", 
               feature2 = "pct_ribo", 
               group.by = "DropletQC") +
  ggtitle("WT - Unfiltered")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_s,
               feature1 = "nFeature_RNA", 
               feature2 = "nCount_RNA", 
               group.by = "DropletQC") +
  ggtitle("WT - Unfiltered")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_s,
               feature1 = "nFeature_RNA", 
               feature2 = "pct_mito", 
               group.by = "DropletQC") +
  ggtitle("WT - Unfiltered")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_s,
               feature1 = "nFeature_RNA", 
               feature2 = "pct_ribo", 
               group.by = "DropletQC") +
  ggtitle("WT - Unfiltered")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_s,
               feature1 = "nuclear_fraction", 
               feature2 = "nCount_RNA", 
               group.by = "DropletQC") +
  ggtitle("WT - Unfiltered")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_s,
               feature1 = "nuclear_fraction", 
               feature2 = "nFeature_RNA", 
               group.by = "DropletQC") +
  ggtitle("WT - Unfiltered")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_s,
               feature1 = "nuclear_fraction", 
               feature2 = "pct_mito", 
               group.by = "DropletQC") +
  ggtitle("WT - Unfiltered")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_s,
               feature1 = "nuclear_fraction", 
               feature2 = "pct_ribo", 
               group.by = "DropletQC") +
  ggtitle("WT - Unfiltered")
dev.off()

################################################################################

# Violin plots
pdf(file.path(WT_plots, "WT_unfiltered_vln_plots_1.pdf"), height = 6, width = 8)
VlnPlot(WT_s, 
        features = "nCount_RNA", 
        group.by = "orig.ident") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # 
VlnPlot(WT_s, 
        features = "nFeature_RNA", 
        group.by = "orig.ident") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # 
VlnPlot(WT_s, 
        features = "pct_mito", 
        group.by = "orig.ident") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # 
VlnPlot(WT_s, 
        features = "pct_ribo", 
        group.by = "orig.ident") + 
  theme(legend.position="none")
dev.off()

# Violin plots grouped by DoubletQC classifications
pdf(file.path(WT_plots, "WT_unfiltered_vln_plots_2.pdf"), height = 6, width = 12)
VlnPlot(WT_s, 
        features = "nCount_RNA", 
        group.by = "DropletQC") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_s, 
        features = "nFeature_RNA", 
        group.by = "DropletQC") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_s, 
        features = "pct_mito", 
        group.by = "DropletQC") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_s, 
        features = "pct_ribo", 
        group.by = "DropletQC") + 
  theme(legend.position="none")
dev.off()


pdf(file.path(WT_plots, "WT_unfiltered_vln_plots_3.pdf"), height = 8, width = 15)
VlnPlot(WT_s, 
        features = "nCount_RNA", 
        group.by = "WT_SC_1") + 
  theme(legend.position="none")
# # # # # # # # # # # # # #
VlnPlot(WT_s, 
        features = "nFeature_RNA", 
        group.by = "WT_SC_1") + 
  theme(legend.position="none")
# # # # # # # # # # # # # #
VlnPlot(WT_s, 
        features = "pct_mito", 
        group.by = "WT_SC_1") + 
  theme(legend.position="none")
# # # # # # # # # # # # # #
VlnPlot(WT_s, 
        features = "pct_ribo", 
        group.by = "WT_SC_1") + 
  theme(legend.position="none")
# # # # # # # # # # # # # #
VlnPlot(WT_s, 
        features = "nuclear_fraction", 
        group.by = "WT_SC_1") + 
  theme(legend.position="none")
dev.off()


# Combine violin plot pdfs
pdf_combine(input = c(file.path(WT_plots, "WT_unfiltered_vln_plots_1.pdf"),
                      file.path(WT_plots, "WT_unfiltered_vln_plots_2.pdf"),
                      file.path(WT_plots, "WT_unfiltered_vln_plots_3.pdf")),
              output = file.path(WT_plots, "0_WT_unfiltered_vln_plots.pdf"))

file.remove(c(file.path(WT_plots, "WT_unfiltered_vln_plots_1.pdf"),
              file.path(WT_plots, "WT_unfiltered_vln_plots_2.pdf"),
              file.path(WT_plots, "WT_unfiltered_vln_plots_3.pdf")))

################################################################################

# UMAP plots
pdf(file.path(WT_plots, "0_WT_unfiltered_UMAPs.pdf"), height = 7, width = 12)
UMAPPlot(WT_s, group.by = "WT_SC_1", label =TRUE) + ggtitle("WT Unfiltered | Clusters")
UMAPPlot(WT_s, group.by = "DropletQC") + ggtitle("WT Unfiltered | DropletQC Classifications")
FeaturePlot(WT_s, features = "nCount_RNA") + ggtitle("WT Unfiltered | # of Counts")
FeaturePlot(WT_s, features = "nFeature_RNA") + ggtitle("WT Unfiltered | # of Genes")
FeaturePlot(WT_s, features = "pct_mito") + ggtitle("WT Unfiltered | % of Mito Genes")
FeaturePlot(WT_s, features = "pct_ribo") + ggtitle("WT Unfiltered | % of Ribo Genes")
FeaturePlot(WT_s, features = "nuclear_fraction") + ggtitle("WT Unfiltered | Nuclear Fraction")
dev.off()
get.all.markers <- function (seurat, ident, path, fname) {
  
  DefaultAssay(seurat) <- "RNA"
  seurat <- ScaleData(seurat, features = row.names(seurat))
  
  Idents(seurat) <- ident
  
  all_markers <- FindAllMarkers(seurat, logfc.threshold = 0.1)
  
  if (is.factor(seurat@meta.data[,ident])) {
    clusters <- levels(seurat@meta.data[,ident])
  } else if (is.character(seurat@meta.data[,ident])) {
    clusters <- sort(unique(seurat@meta.data[,ident]))
  }
  
  all_markers_list <- vector(mode="list", length = length(clusters))
  names(all_markers_list) <- clusters
  
  for (clust in clusters) {
    all_markers_list[[clust]] <- all_markers[which(all_markers$cluster == clust), 
                                             c("gene", "p_val", "avg_log2FC", 
                                               "pct.1", "pct.2", "p_val_adj")]
  }
  write.xlsx(all_markers_list, 
             file.path(path, fname), 
             rowNames = FALSE)
  
  return(all_markers_list)
  
}
all_markers <- get.all.markers(seurat = WT_s,
                               ident = "seurat_clusters",
                               path = WT_markers,
                               fname = "0_WT_cluster_markers.xlsx")
########################### WT Low-quality Filtering ###########################

WT_f <- subset(WT_s, 
               subset = nFeature_RNA <= 8000 &
                 nFeature_RNA >= 1500 & 
                 nCount_RNA >= 1000 &
                 nCount_RNA <= 40000 &
                 pct_mito <= 25)

Idents(WT_f) <- "DropletQC"
WT_f <- subset(WT_f, idents = "cell")

# Number of cells
nrow(WT_f@meta.data) # 5524
########################### WT Doublet Classifications #########################

WT_f <- FindVariableFeatures(WT_f, 
                             selection.method = "vst", 
                             nfeatures = 2000)
WT_f <- ScaleData(WT_f)
WT_f <- RunPCA(WT_f)

# Determine number of dimensions that capture at least 80% of the variation
pct <- WT_f[["pca"]]@stdev / sum(WT_f[["pca"]]@stdev) * 100
cumu <- cumsum(pct)
PC_80_pct <- which(cumu >= 80)[1]
PC_80_pct #33

# UMAP and clustering
WT_f <- FindNeighbors(WT_f, 
                      dims = 1:PC_80_pct)
WT_f <- RunUMAP(WT_f, 
                dims=1:PC_80_pct)

WT_f <- FindClusters(WT_f, 
                     resolution = 0.6)

WT_f <- AddMetaData(WT_f, WT_f$seurat_clusters, col.name = "WT_SC_2")

################################################################################

# Doublet detection using scDblFinder
WT_sce <- as.SingleCellExperiment(WT_f, 
                                   assay = "RNA")

set.seed(123)

WT_sce <- scDblFinder(WT_sce, 
                      clusters="WT_SC_2")

identical(WT_sce@colData@rownames, row.names(WT_f@meta.data)) # TRUE
identical(WT_sce$WT_SC_2, WT_f@meta.data[,"WT_SC_2"]) # TRUE

WT_f <- AddMetaData(WT_f, WT_sce$scDblFinder.class, col.name = "scDblFinder.class")
WT_f <- AddMetaData(WT_f, WT_sce$scDblFinder.score, col.name = "scDblFinder.score")

saveRDS(WT_f, file.path(WT_rds, "1_WT_seurat.rds"))
###################### WT Plots w/ scDblFinder Classifications #################

pdf(file.path(WT_plots, "1_WT_doublet_scatter_plots.pdf"), height = 8, width = 12)
FeatureScatter(WT_f,
               feature1 = "nCount_RNA", 
               feature2 = "nFeature_RNA", 
               group.by = "scDblFinder.class")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_f,
               feature1 = "nCount_RNA", 
               feature2 = "pct_mito", 
               group.by = "scDblFinder.class")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_f,
               feature1 = "nCount_RNA", 
               feature2 = "pct_ribo", 
               group.by = "scDblFinder.class")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_f,
               feature1 = "nFeature_RNA", 
               feature2 = "nCount_RNA", 
               group.by = "scDblFinder.class")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_f,
               feature1 = "nFeature_RNA", 
               feature2 = "pct_mito", 
               group.by = "scDblFinder.class")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_f,
               feature1 = "nFeature_RNA", 
               feature2 = "pct_ribo", 
               group.by = "scDblFinder.class")
dev.off()

################################################################################

# Violin plots grouped by scDblFinder classifications
pdf(file.path(WT_plots, "1_WT_doublet_vln_plots.pdf"), height = 6, width = 12)
VlnPlot(WT_f, 
        features = "nCount_RNA", 
        group.by = "scDblFinder.class") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "nFeature_RNA", 
        group.by = "scDblFinder.class") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "pct_mito", 
        group.by = "scDblFinder.class") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "pct_ribo", 
        group.by = "scDblFinder.class") + 
  theme(legend.position="none")
dev.off()

################################################################################

# Violin plots grouped by scDblFinder clusters
pdf(file.path(WT_plots, "1_WT_vln_scdblfinder_clusters.pdf"), height = 6, width = 12)
VlnPlot(WT_f, 
        features = "nCount_RNA", 
        group.by = "WT_SC_2") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "nFeature_RNA", 
        group.by = "WT_SC_2") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "pct_mito", 
        group.by = "WT_SC_2") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "pct_ribo", 
        group.by = "WT_SC_2") + 
  theme(legend.position="none")
dev.off()

################################################################################

# Violin plots grouped by original clusters
pdf(file.path(WT_plots, "1_WT_vln_original_clusters.pdf"), height = 6, width = 12)
VlnPlot(WT_f, 
        features = "nCount_RNA", 
        group.by = "WT_SC_1") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "nFeature_RNA", 
        group.by = "WT_SC_1") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "pct_mito", 
        group.by = "WT_SC_1") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "pct_ribo", 
        group.by = "WT_SC_1") + 
  theme(legend.position="none")
dev.off()

################################################################################

# UMAP plots
pdf(file.path(WT_plots, "1_WT_doublet_UMAPs.pdf"), height = 7, width = 12)
UMAPPlot(WT_f, group.by = "WT_SC_2", label =TRUE) + ggtitle("WT Doublets | Clusters")
UMAPPlot(WT_f, group.by = "scDblFinder.class") + ggtitle("WT Doublets | DropletQC Classifications")
FeaturePlot(WT_f, features = "nCount_RNA") + ggtitle("WT Doublets | # of Counts")
FeaturePlot(WT_f, features = "nFeature_RNA") + ggtitle("WT Doublets | # of Genes")
FeaturePlot(WT_f, features = "pct_mito") + ggtitle("WT Doublets | % of Mito Genes")
FeaturePlot(WT_f, features = "pct_ribo") + ggtitle("WT Doublets | % of Ribo Genes")
FeaturePlot(WT_f, features = "scDblFinder.score") + ggtitle("WT Doublets | scDblFinder Score")
dev.off()
############################### WT Doublet Removal #############################

# Filter out the cells classified as doublets
Idents(WT_f) <- "scDblFinder.class"
WT_f <- subset(WT_f, idents = "singlet")

# Number of cells
nrow(WT_f@meta.data) # 5055

################################################################################

WT_f <- FindVariableFeatures(WT_f, 
                             selection.method = "vst", 
                             nfeatures = 2000)
WT_f <- ScaleData(WT_f)
WT_f <- RunPCA(WT_f)

# Determine number of dimensions that capture at least 80% of the variation
pct <- WT_f[["pca"]]@stdev / sum(WT_f[["pca"]]@stdev) * 100
cumu <- cumsum(pct)
PC_80_pct <- which(cumu >= 80)[1]
PC_80_pct #33

# UMAP and clustering
WT_f <- FindNeighbors(WT_f, 
                      dims = 1:PC_80_pct)
WT_f <- RunUMAP(WT_f, 
                dims=1:PC_80_pct)

WT_f <- FindClusters(WT_f, 
                     resolution = 0.6)

WT_f <- AddMetaData(WT_f, WT_f$seurat_clusters, col.name = "WT_SC_3")

saveRDS(WT_f, file.path(WT_rds, "2_WT_seurat.rds"))
############################# WT Cell Cycle Scoring ############################

s.genes <- c("Mcm5","Pcna","Tyms","Fen1","Mcm2","Mcm4","Rrm1","Ung","Gins2",
             "Mcm6","Cdca7","Dtl","Prim1","Uhrf1","Hells","Rfc2","Rpa2","Nasp",
             "Rad51ap1","Gmnn","Wdr76","Slbp","Ccne2","Ubr7","Pold3","Msh2",
             "Atad2","Rad51","Rrm2","Cdc45","Cdc6","Exo1","Tipin","Dscc1","Blm",
             "Casp8ap2","Usp1","Clspn","Pola1","Chaf1b","Brip1","E2f8")

g2m.genes <- c("Hmgb2","Cdk1","Nusap1","Ube2c","Birc5","Tpx2","Top2a","Ndc80",
               "Cks2","Nuf2","Cks1b","Mki67","Tmpo","Cenpf","Tacc3","Pimreg",
               "Smc4","Ccnb2","Ckap2l","Ckap2","Aurkb","Bub1","Kif11","Anp32e",
               "Tubb4b","Gtse1","Kif20b","Hjurp","Cdca3","Jpt1","Cdc20","Ttk",
               "Cdc25c","Kif2c","Rangap1","Ncapd2","Dlgap5","Cdca2","Cdca8",
               "Ect2","Kif23","Hmmr","Aurka","Psrc1","Anln","Lbr","Ckap5",
               "Cenpe","Ctcf","Nek2","G2e3","Gas2l3","Cbx5","Cenpa")

WT_f <- CellCycleScoring(WT_f,
                         s.features = s.genes,
                         g2m.features = g2m.genes,
                         set.ident = FALSE)

saveRDS(WT_f, file.path(WT_rds, "3_WT_seurat.rds"))
############################## WT Filtered Plots ###############################

pdf(file.path(WT_plots, "2_WT_scatter_plots.pdf"), height = 8, width = 12)
FeatureScatter(WT_f,
               feature1 = "nCount_RNA", 
               feature2 = "nFeature_RNA", 
               group.by = "orig.ident")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_f,
               feature1 = "nCount_RNA", 
               feature2 = "pct_mito", 
               group.by = "orig.ident")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_f,
               feature1 = "nCount_RNA", 
               feature2 = "pct_ribo", 
               group.by = "orig.ident")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_f,
               feature1 = "nFeature_RNA", 
               feature2 = "nCount_RNA", 
               group.by = "orig.ident")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_f,
               feature1 = "nFeature_RNA", 
               feature2 = "pct_mito", 
               group.by = "orig.ident")
# # # # # # # # # # # # # # # # # # # # # # # #
FeatureScatter(WT_f,
               feature1 = "nFeature_RNA", 
               feature2 = "pct_ribo", 
               group.by = "orig.ident")
dev.off()

################################################################################

# Violin plots grouped by scDblFinder classifications
pdf(file.path(WT_plots, "2_WT_vln_plots.pdf"), height = 6, width = 12)
VlnPlot(WT_f, 
        features = "nCount_RNA", 
        group.by = "orig.ident") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "nFeature_RNA", 
        group.by = "orig.ident") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "pct_mito", 
        group.by = "orig.ident") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "pct_ribo", 
        group.by = "orig.ident") + 
  theme(legend.position="none")
dev.off()

################################################################################

# Violin plots grouped by final clusters
pdf(file.path(WT_plots, "2_WT_vln_final_clusters.pdf"), height = 6, width = 12)
VlnPlot(WT_f, 
        features = "nCount_RNA", 
        group.by = "WT_SC_3") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "nFeature_RNA", 
        group.by = "WT_SC_3") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "pct_mito", 
        group.by = "WT_SC_3") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "pct_ribo", 
        group.by = "WT_SC_3") + 
  theme(legend.position="none")
dev.off()

################################################################################

# Violin plots grouped by original clusters
pdf(file.path(WT_plots, "2_WT_vln_original_clusters.pdf"), height = 6, width = 12)
VlnPlot(WT_f, 
        features = "nCount_RNA", 
        group.by = "WT_SC_1") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "nFeature_RNA", 
        group.by = "WT_SC_1") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "pct_mito", 
        group.by = "WT_SC_1") + 
  theme(legend.position="none")
# # # # # # # # # # # # # # # # # # # # # # # #
VlnPlot(WT_f, 
        features = "pct_ribo", 
        group.by = "WT_SC_1") + 
  theme(legend.position="none")
dev.off()

################################################################################

# UMAP plots
pdf(file.path(WT_plots, "2_WT_UMAPs.pdf"), height = 7, width = 12)
UMAPPlot(WT_f, group.by = "WT_SC_3", label =TRUE) + ggtitle("WT | Clusters")
UMAPPlot(WT_f, group.by = "Phase", label =FALSE) + ggtitle("WT | Cell Cycle")
FeaturePlot(WT_f, features = "nCount_RNA") + ggtitle("WT | # of Counts")
FeaturePlot(WT_f, features = "nFeature_RNA") + ggtitle("WT | # of Genes")
FeaturePlot(WT_f, features = "pct_mito") + ggtitle("WT | % of Mito Genes")
FeaturePlot(WT_f, features = "pct_ribo") + ggtitle("WT | % of Ribo Genes")
dev.off()
WT_f$RNA_snn_res.0.6 <- NULL 
WT_f$seurat_clusters <- NULL 

saveRDS(WT_f, file.path(getwd(), "WT_seurat.rds"))

Page 1 Page 2 Page 3 Page 4 Page 5 Page 6 Page 7